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Abstract 
> 

(^ . The Vlasov-Nordheim equation is solved numerically on a lattice for 

t>. 

nuclear matter in two dimensions. We discuss the reliability of the model 

at normal density and then study the response of the system to small 

perturbations. We find deterministic chaos inside the spinodal zone 

O, 

where fragment formation occurs. We discuss in detail the dynamical 

features of this phenomenon in order to clarify the mechanisms leading 

to nuclear disassembly in heavy-ion collisions. 
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Many recent experimentsa in heavy-ion collisions at energies above 50 MeV/A 
have shown a very rapid break-up of the hot composite system formed in the reac- 
tion into several big fragments with Z > 3. This disassembly of the nuclear system, 
which can include a substantial fraction of the colliding masses, falls under the generic 
name of " nuclear multifragmentation" . In parallel with experiments there has been 
an intense theoretical effort in the understanding of the mechanisms underlying the 
phenomenonH. In fact a detailed and systematic study is expected to provide fun- 
damental properties of nuclear matter. However, though many models are able to 
reproduce the data, the origin of multifragmentation has not been sufficiently clari- 
fied up to nowa. 

It seems nowadays established that multifragmentation cannot be explained by 
the same scenarios valid at lower energies. The characteristic time ranges from one 
hundred to few hundreds fm/cu. Therefore it is a very fast process in comparison with 
phenomena like fission or compound nucleus formation whose typical time is of the 
order of thousands fm/c. Mechanical instabilities-!! play a fundamental role and they 
should be taken into due account when the excited nuclear system enters the spinodal 
zone of nuclear matter Equation of State (EOS). The latter is the region where the 
compressibility, i.e. the derivative of the pressure with respect to the density, becomes 
negative and the system is unstable to small perturbations. 

The undoubted success of the statistical modelsnIJ, in explaining some of the mul- 
tifragmentation features, indicates that the phase space dominates the population of 
the different final channels. In such a fast process, the tendency of filling uniformly the 
phase space cannot occur in each collisional event. The system has not enough time 
to explore, during the reaction, the whole phase space, as can happen for compound 
nucleus collisions. It can be expected, therefore, that the reaction dynamics is dom- 
inated by the phase space only when the physical quantities are averaged over large 



sets of events. This assumption implies that the nuclear dynamics in the multifrag- 
mentation regime is irregular or chaotic enough to produce, at least approximately, 
a uniform "a priori" probability to populate each region of the available phase space. 
In particular, the formation of the final fragments must follow an irregular dynamics. 
This conjecture can be also inferred by the large event-by-event fluctuations observed 
experimentally on the charge and mass distributions^. 

Within this scenario, however, it is still not clear at which stage the dynamics is 
dominated by a chaotic behaviour and which is the mechanism able to produce the 
expected phase space mixing. 

In the present work the problem of fragment formation is approached in a sim- 
plified way. It is assumed that the nuclear system after an initial compression slowly 
expands entering into the spinodal region. Some evidences of such a stage of the re- 
action in central heavy ion collisions have been found in computer simulations based 
on the BUU schemeQ'B. The fragments are supposed to emerge from this expanded 
nuclear system through the unstable density fluctuations of the spinodal zone. We 
study in details this phase of the nuclear dynamics and show that it exhibits a chaotic 
behavior. In section 1 we discuss a simple example to illustrate the main features 
of deterministic chaos. The theoretical background for the nuclear case is based on 
the Vlasov equation, which is overviewed in section 2. In section 3 we discuss the 
numerical method adopted to solve the Vlasov equation on a lattice and study the 
corresponding two-dimensional Equation of State (EOS). The mean field dynamics is 
investigated in section 4. The density profile of the initially almost uniform nuclear 
system is followed until fragments can be identified. The dynamics of this process is 
analyzed in order to investigate the possible non-linear and chaotic behaviour, with 
methods analogous to the ones well known in the theory of dynamical systems'! To 
this purpose, a distance between "trajectories" is introduced, and the analysis in 



terms of Lyapunov exponents is performed. In section 5 we discuss the consequences 
of a non-linear and erratic dynamics for nuclear multifragmentation. Finally in sec- 
tion 6 we draw our conclusions. The results reported in this paper are a complete 
and detailed review of an investigation started along these lines in refs.QH. 

1. Deterministic chaos in a simple example. 

In the last two decades the study of low-dimensional dynamical systems has given 
a fundamental contribution to the understanding of the onset of irregular motion. Due 
to the nonintegrability of the sets of deterministic equations an erratic and unpre- 
dictable behaviour, deterministic chaos, can follow. What was previously considered 
as spurious noise has found rigorous mathematical and experimental foundations. 
This revolutionary concept has changed drastically our view of classical mechanics 
and is revealing important implications also in quantum physics. For a general review 
see ref.H The consequences of this new paradigm which is shading light into apparently 
unrelated different fields are difficult to overcast and will probably be an important 
guideline of research for the next decades. In this work, we explore the dynamics of 
nuclear matter from this new perspective. But before going to the discussion of the 
main topic, we will consider the dynamics of a very simple system. This example is 
used to illustrate the main features of a chaotic unstable system and the conditions 
under which deterministic chaos can occur. 

Let us consider the two-dimensional problem described by the hamiltonian 

n = \{v\ + vl) + ql + ql - ^q\q\ (i) 

which corresponds to two coupled harmonic oscillators with unit mass. It can be easily 
seen that, for a given total energy E, the available space has a fourfold structure, with 
four "branches" extending to infinity. If \x ^ 0, the dynamics of the particle inside the 



potential is chaotic. In other words, it is impossible to find an analytical solution for 
the equations of motion corresponding to the Hamiltonian (1). We have two degrees 
of freedom and only one constant of the motion, i.e. the total energy. This implies 
an extreme sensitivity of the dynamics on the initial conditions and a very irregular 
motion. 

One can consider a particle starting its motion inside this potential or impinging 
on it from outside. In the latter case the chaoticity of the interaction region manifests 
itself in strong irregularities, at all scales, of the final observables - final scattering 
angle, final internal excitation, etc. - as a function of the initial conditions. We have 
the well known phenomenon of chaotic scattering^. An incoming particle can be 
trapped for a long period inside the interaction region, where it bounces erratically 
to and fro before escaping. The extreme internal sensitivity influences drastically the 
final result. 

Classical chaotic scattering is therefore intimately connected with the presence of 
trajectories which remain trapped in the interaction region for a time long enough 
that chaos can set in. 

These considerations can be readily extended to the case of trajectories which 
starts directly from the interaction region and which therefore will describe systems 
which are unstable and decay. In particular, the Hamiltonian (1) does not describe 
a proper scattering situation, since the potential diverges at infinity. However it can 
describe a system which escapes from the interaction region after some trapping time. 
Since the system is chaotic, again the trapping time and the direction of emission 
depend in a very irregular way on the initial conditions. This point is illustrated in 
fig.l. For each trajectory starting from the internal region, one of the coordinates of 
the particle is recorded. The distance at which we take this value is so large to be 
considered asymptotic. The particle cannot return anymore to the interaction region. 



In the plot, the value of the final coordinate qj is reported as a function of the initial 
one qi. It is important to notice that the considered set of initial conditions lies on the 
same energy surface. The observed irregular structure is very similar to the one usual 
encountered for the deflection function in chaotic scattering systems. If chaoticity is 
strong enough, the distribution of the trapping time is actually exponential^. More 
formally, if we call T the average trapping time and r = 1/A the average rate of 
divergence between trajectories, being A the average Lyapunov exponent, one can 
write, according to TelEJ 

T = r/{l-D) (2) 

where D is the average linear Hausdorff dimension of the so-called strange repelloirB. 
If D > 0, which is the condition for a genuine chaotic scattering, then T > t. 

This condition can be used as a criterium to argue if a system presents a chaotic 
dynamics or not. It is useful in applications, because it is usually numerically difficult 
to assess the existence of positive D by an extensive direct sampling of the whole 
phase space. Still it must be used with caution, due to the unavoidable uncertainties 
in the numerical values of the Lyapunov exponents, and therefore it can give only an 
indication of a possible chaotic dynamics. 

The analysis can be extended to models with many degrees of freedom. In par- 
ticular, for a fluid system the relevant degrees of freedom can be identified with its 
linear eigenmodes, and in the harmonic approximation a fluid is exactly equivalent 
to a set of harmonic oscillators. Such a system, which corresponds to the case /i = 
in eq. (1), is integrable, since the energy of each mode is conserved. Therefore, in 
this limit a fluid can undergo only a regular dynamical evolution. If the amplitude of 
the perturbation is large enough, the modes of the fluid can be coupled between each 
other, as in eq. (1) for fi ^ 0. In this case the problem is, in general, non-integrable 
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and chaotic dynamics can be present. 

If a uniform homogeneous fluid is initially in its spinodal region, where the com- 
pressibility is negative, it will spontaneously "escape" towards a non-homogenous 
phase, where droplets are formed. However, before the final droplets are formed, the 
fluid can spend enough time in the spinodal region to allow non-linearity and chaotic 
motion to set in and dominate the time evolution. The situation is quite analogous 
to the example of eq. (1), where the particle escapes outside the interaction region 
but has enough time to experience the chaotic dynamics present in the inner sector 
of the potential, as illustrated in fig.l. If chaoticity is really present, then the av- 
erage Lyapunov exponent A must be positive and, as discussed above, the criterium 
T > t should be satisfied. Here T is the characteristic time for droplet formation 
and t = 1/A is the average time of divergence between "trajectories" of the fluid. In 
the next sections we study - within the semiclassical approximation - the dynamics of 
nuclear matter inside the spinodal region. We will show, following a reasoning simi- 
lar to the one discussed above, that nonlinearities and positive Lyapunov exponents 
are found inside the spinodal region. Therefore deterministic chaos characterizes the 
formation of fragments. 

2. Theoretical background 

2.1. The Vlasov equation 

In the multifragmentation energy range the dynamics is characterized by an inter- 
play between a purely mean field evolution, typical of a low-energy phenomenology 
(fusion, inelastic reactions), and two-body collisions due to the partial relaxation 
of the Pauli principle. A possible way of incorporating both of these effects is the 
inclusion of the residual interaction in models like TDHF, but this procedure is diffi- 
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cult to carry on numerically in realistic calculations!.! For this purpose semiclassical 
methods have been developed, whose starting point is the Wigner transform of the 
time-dependent Hartree-Fock (TDHF) equation for the one-body density matrix. By 
neglecting powers of h higher than two^, one obtains the usual Vlasov equation which 
reads 

% + ~^rf-^rU[p(r)}.V p f = (3) 

/(r, p; t) is the Wigner phase-space distribution function, (r, p) are the space and 
momentum coordinates and U is a self-consistent single-particle potential depending 
on the density p = J dp f . When eq. (3) includes a Nordheim-type collision integral 
on the right-hand side, I[f], it is generally referred in literature as BUU (Boltzmann- 
Uehling-Uhlenbeck) or VUU (Vlasov-Uehling-Uhlenbeck) or LV (Landau- Vlasov) or 
VN (Vlasov-Nordheim) equation. For a complete review, see ref.li-3. The collision 
integral reads 

^[/](r,Pi,£) = / dp2 dpi/ dp 2 > Ui'h'hh - f&fv f2')vvii'y (4) 

where fj = f(sj,t) is the phase-space occupation probability at the location Sj = 
(r,pj), and fj — 1 — fj is the Pauli blocking factor of the final states. iui 2 i'2' is the 
microscopic transition rate for the collision vertex (P1P2) — ► (pvPw), and is related 
to the nucleon-nucleon scattering cross section. This collision term describes the 
momentum changes of two interacting particles during a collision with blocking factors 
forbidding transitions leading to occupied final states. 

We notice the formal identity of eq. (3) with the Liouville equation of classical me- 
chanics for a fluid of non-interacting particles. This underlines its intrinsic classical 
character and guarantees that quantum effects like the Pauli principle are conserved. 
It can be easily demonstrated that the Vlasov-Nordheim equation satisfies the conser- 
vation laws of mass, momentum and energy. The presence of an average potential U 



dependent on the density introduces a non-linear constraint on the solution of eq.(3). 
This is a common point to all mean field theories. 

Because of that, the Vlasov-Nordheim equation is difficult to solve numerically, 
in spite of its apparent simplicity. A large amount of literature has been devoted to 
its numerical solution techniques in the case of systems interacting through Coulomb 
fields like plasmas, to which the Vlasov equation was first applied. 

Later on we will discuss the implications of non-linearity on physical processes 
like multifragmentation. 

2.2. Linear response for the Vlasov equation 

The linearized version of the Vlasov equation has been successfully applied to the 
study of small amplitude oscillations, see ref.E-3't3, and it represents the starting point 
for a phase space approach to RPA. For small variations of the distribution function 
around the equilibrium solution f (r,p 2 ) 

/(r,p,£) = / (r,p 2 )+2(r,p,i) (5) 

we can expand eq.(3) to the first order in g. A Thomas- Fermi approximation for the 
ground state - which is obviously a solution of eq. (3) - can be used. That is 

*'<■'•& = <&<f°&-&> (6) 

where D is the dimension of the physical space. Then one gets an equation for 
the Fourier transform of g(r,p,t) which can be analytically solved and produces a 
dispersion relation for the frequency ui for each wave number k. For stable systems 
the linear response theory yields real values for the energies u = u(k) of the normal 
modes, which can be represented by plane waves propagating in opposite directions. 
When the system is unstable, some energies are complex and the imaginary part 
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is related to the time growth of instabilities. In both cases, the dependence of the 
energy u on k is linear in this approximation, with a correction factor depending on 
the equation of state and the density of the system, see ref.l^fi3'E-l We like to remind 
that the linear response approximation has a limited range of validity, since it can 
be applied only to small amplitude variations of the initial reference state. Therefore 
its applicability to processes like multifragmentation, where large fluctuations are 
involved, might be not completely correct, as it will be discussed later. 

3. The model 

3.1. The lattice calculation 

As previously explained in the introduction, we study the behavior of a dilute 
nuclear system whose density fluctuations are growing inside the spinodal region. In 
order to get numerically robust results, we solved the Vlasov-Nordheim equation on a 
lattice, using the same code of refill, but neglecting the stochastic contribution to the 
collision integral. We divide the single particle phase space into several small cells, 
each of volume V D = Ar D ■ Ap D , being Ar and Ap the cell sidelengths respectively in 
coordinate and momentum space. Most numerical problems in calculating the Vlasov 
evolution of nuclear systems arise from the need to smooth the one-body density. For 
this purpose we must employ a lattice with a big number of small cells in order to get a 
nice paving of the phase space. Typically Ar is 1 fm or less, while Ap should be smaller 
than the Fermi momentum. Of course the discretization introduces some numerical 
error on the physical observables, but it gets smaller and smaller as we decrease the 
cell size and increase their number. Therefore the main limit of the lattice method is 
the memory resources of the computers and the huge computing time requested. For 
this reason, we performed all the calculations on a two-dimensional lattice. 
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We have studied a fermion gas situated on a large torus with periodic boundary 
conditions, and its size is kept constant during the evolution. The torus sidelengths 
are equal to L x = 51 fm and L y = 15 fm. We employed in momentum space 51x51 
small cells of size Ap x = Ap y = 40 MeV/c, while in coordinate space Ax = 0.3333 fm 
and Ay = 15 fm, i.e. we have only one big cell on the y-direction. 

The initial local momentum distribution was assumed to be the one of a Fermi 
gas at a fixed temperature T. We employ a local Skyrme interaction which generates 
a mean field U[p] = t (p/po) + £3 (p/po) 2 . The density p is folded along the re- 
direction with a gaussian of width p = 0.61/m, in order to give a finite range to 
the interaction. The parameters of the force to an d t 3 have been chosen in order to 
reproduce correctly the binding energy of nuclear matter at zero temperature, and 
this gives t = —100.3 MeV and t 3 = 48 MeV. The resulting EOS gives a saturation 
density in two dimensions equal to po = 0.55 fm" 2 which corresponds to the usual 
three-dimensional Fermi momentum equal to Pp = 260 MeV jc. 

Then a complete dynamical evolution is performed by subdividing the total time 
in small time steps, each equal to At = 0.5 fm/c. 

For more details concerning the mean field propagation on the lattice and the 
exact calculation of the collision integral, the reader is referred to refill 

3.2. Two-dimensional nuclear matter equation of state 

We calculate the nuclear matter equation of state (EOS) for an homogeneous 
two dimensional system. We mainly follow the definitions given in ref.H3 for the 
three dimensional case. For a gas of particles interacting through a local Skyrme 
single-particle potential, U[p], we calculate the density of the free energy F and the 
corresponding pressure P. They read 

F = H-TS (7) 
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P = -F + p(p + U\p]) (8) 

In eq.(7) His the spatial energy density, and it is given by summing the kinetic energy 
density K and the potential energy density E — J U[p]dp. p denotes the chemical 
potential, T is the temperature and S is the density of entropy, for which we use the 
standard definition for noninteracting particlestJ. 
If we define the parameter 

V = % (9) 

the single-particle density, the kinetic energy density and the entropy density can be 
expressed through the Fermi integrals Jo and J\ 

P=—,Jo(v) (10) 

Tin 

K=—,J 1 ( V ) (11) 

Tin 

Am 

S = —5-JM-VP (12) 

Tin 1 

where the Fermi integral is given by 

r e u de 
JM = / t^^ (13) 

J 1 + er 'I 

Since the integral J$ is analytical, the chemical potential p is readily calculated from 
the single-particle density p 

p = T log[exp{^£) - 1] (14) 
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From the above definitions, the expressions for the density of the free energy and the 
pressure are 

2T72 

F = ~— 2 Ji(v)+E + fxp (15) 

Tin 

2T72 

P =^tf JM-E + pU\p] (16) 

In Fig. 2 we show the two-dimensional nuclear matter equation of state (EOS), cal- 
culated for the Skyrme interaction we discussed in the previous subsection. In the 
part a) of the figure, the free energy per particle is displayed along the isotherms for 
different values of the temperature T as function of the density, while in part b) the 
corresponding thermostatic pressure is shown. This recalls the qualitative features of 
a classical Van der Waals gas. We notice that the general trend of a two-dimensional 
EOS does not differ appreciably from the realistic three-dimensional casell2l. 

It has to be noticed that the pressure P can be also calculated by the thermo- 
dynamical relationship P = —p 2 ^-, being F the free energy per particle. It can be 
readily verified that this procedure gives the same expression of eq.(16). The latter 
result indicates that in our particular simple model of nuclear matter, the thermody- 
namical relationships are exactly satisfied. This is a consequence of the simple form 
adopted for the interaction. 

The spinodal zone can be identified in fig. 2b) as the region where the compress- 
ibility , i.e. the slope of the pressure P as a function of the density p, is negative. At 
increasing temperature the density interval where |p is negative becomes smaller and 
smaller. Finally at the critical temperature T c for the liquid-gas phase transition the 
spinodal zone reduces to a point where the corresponding isotherm has an inflection. 
For T > T c the pressure P is a monotonic increasing function of p. From the figure 
one can deduce T r ~ 16 MeV. 



13 



4. Mean field dynamics 

In this section we want to study in detail the response of nuclear matter to small 
perturbations. We neglect for the moment the collision term 1(f) and solve numerically 
only the Vlasov equation. If non-linear terms are negligible, linear response theory 
gives an accurate description of the time evolution. However we will see that this is 
not always the case. 

Let us consider an initial sinusoidal perturbation along the x— direction in the 
density profile. That is 

p{x) = [ 1 + Ssin(kx)] p , (17) 

with 

k=-jr±, (18) 

being L x the size of the torus along the x-direction and p = J dxp(x)/L x . We study 
the time evolution of this profile considering a small amplitude for the perturbation, 
i.e. 5 = 0.01. In fig.3 and 4 we show this evolution for p/po — 0.8 for the cases 
rik = 4, 5. In these figures the mode initially excited is damped in a time range which 
is around 40 fm/c. This is the usual Landau damping of Fermi liquid theory which is 
a well known linear phenomenon. Actually, since we are working in a discrete lattice, 
we have also a small initial excitation of a subharmonic of higher order, i.e. nk = 12 
and rik = 25, which is also damped very rapidly. In a previous paperQ we had shown 
a similar damping but at normal density. In that case the damping is slower and 
occurs within 70 fm/c. The Landau damping for nuclear matter in two dimensions is 
an important issue by itself and will be discussed in detail in a forthcoming paperEa 
The behavior shown in the two previous figures indicates that initial perturbations 
evolve linearly in time outside the spinodal zone. Moreover a small difference in the 
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initial conditions does not produce a different dynamical evolution, as we will discuss 
later. 

The situation is completely different if one analyses the time evolution inside the 
spinodal zone, i.e. for p/po < 2/3. In fig. 5 we display the evolution of three initial 
perturbations at p/po = 0.4. The small initial modulations evolve and form fragments 
with different shapes and sizes. In our case, due to the fact that we are considering 
a two-dimensional torus of nuclear matter, no real fragmentation occurs and we call 
"fragments" the macroscopic structures of the density profile. The behavior shown 
in fig. 5 is typical of a dynamical system in a chaotic region, where very small initial 
perturbations are rapidly amplified and distorted. The relevant characteristic is that 
the density profiles do not increase rapidly in their amplitude only, they also change 
their shape strongly during the time evolution. These features are strong indications 
of a non-linear regime. 

In fig. 6 we show the power spectra corresponding to the modes involved in fig. 5. 
These spectra are different in each case and illustrate a population of more and more 
modes as the time evolution goes on. Hence the behaviour of the system for p/po = 0.4 
show a sensitive dependence on the initial conditions and an unpredictability of the 
evolution for a given initial perturbation. Small initial differences will produce large 
final deviations^. A similar feature was illustrated in ref.0 for p/po = 0.5. As we will 
clarify in the following this behavior is typical inside the spinodal region and does 
not depend on the particular choice of varying the initial conditions. In the present 
case we have changed the initial nk, but we could have modified either the average 
density or the amplitude of the initial perturbation. 

The sensitivity to the initial conditions occurs, in the dynamical evolution, to- 
gether with a dominant role of the non-linear terms in the Vlasov equation. This is 
proved by the time evolution of the strength of the most important modes as illus- 
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trated in fig. 7. In the latter the numerical calculations are reported as open squares 
in steps of 5 fm/c. Linear fits of the initial evolution are shown for comparison. It 
is evident that the growth is linear only in the first stages of the evolution while it 
strongly deviates from the fits after ~ 25/m/c. In general one cannot trust linear 
response theory after the very first dynamical stages. It is true that some modes 
follow a linear evolution for a longer time, however in general they do correspond to 
very small wavelengths. Therefore they do not affect the size of the main emerging 
fragments but only their surface modulations. One can consider A/2 as the diameter 
of the nuclear fragment formed. In concluding this section one can say that when 
the system lies inside the spinodal region many modes are coupled and they interact 
strongly among each other. The final fragments are the result of this dynamics which 
goes beyond the linearized equations discussed in section 2.2. 



4.1. Calculation of the largest Lyapunov exponent 

The arguments discussed so far provide a strong indication that the dynamics is 
chaotic beyond the spinodal line, but they do not provide a real proof that this is so. 
Therefore in this section we present numerical calculations which define the degree of 
chaoticity according to the usual techniques adopted in the study of low-dimensional 
dynamical systemsO. The discussion follows the same lines of ref.Qu but gives a more 
complete and detailed overview. 

A way to characterize quantitatively the dynamics in a chaotic regime is by means 
of the rate of divergence of two nearby trajectories. The mean rate is usually called 
largest Lyapunov exponents and it is calculated by means of the expression 



A = lim lim \(t) (19) 

t-t-oo d -*0 
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where 

X(t) = log(rf(t)/ ^ o) . (20) 

In the latter d(t) is the distance between two phase space trajectories, along an 
unstable direction. In our notation do = d(0). We have chosen as metric the following 
one 

d(t) = £|p (1 W)-p (2) (^)l/^ c , (21) 

i 

where the index i runs over the N c cells in ordinary space, and p( l >(xi), p( 2 '(xi) are 
the densities in the cell Xi for the trajectories 1 and 2 respectively. The definition of 
eq. (21) represents the difference between the two density profiles and is sufficient 
for our purpose, although possible differences in momentum space are averaged out. 
It should include the contribution of all the unstable modes which dynamically grow 
up during the evolution. 

In fig. 8 we show the time evolution of A(£) for various average densities. Inside 
the spinodal region A(£) converges to a limiting value A, represented by a dashed 
line, which is just the largest Lyapunov exponent. This convergence occurs only for 
a limited time interval, which is the one needed for fragment formation. As it will 
be discussed below this is because we are describing a transient phenomenon. A 
different behavior is observed when the system evolves outside the spinodal region. 
In fact in this case the trajectories do not diverge exponentially and X(t) goes slowly 
towards zero as time increases. This result is analogous to the general case occurring in 
simple dynamical systems, when A finite and positive distinguishes a chaotic behavior 
from a regular one where A = 0. The values of A at a temperature T=3 MeV are 
reported on the second column of table 1. Some of the features here discussed were 
already presented in ref! where actually the evolution was followed for a shorter 
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time. For the calculations presented in this paper we improved the algorithm used 
and we could follow the time evolution for a longer period maintaining an error in the 
energy conservation smaller than 1%. Now one can see that, inside the spinodal zone, 
X(t) remains around the plateau - which indicates the value of the largest Lyapunov 
exponent - only for a limited time range. This fact has a simple explanation. In fact, 
the time at which the plateau stops corresponds to the primary fragment formation. 
Hence after fragments are formed, the different modes interact weaker than before 
among each other, because along the x-direction there are many points where the 
density is almost zero. These fragments stop increasing their size and start to move 
along the x-direction. They eventually collide against each other and therefore, in 
effect, the fragment configuration continues to change. However this dynamical regime 
is different from the previous one and we will not consider the time evolution following 
the plateau. It is important to notice that estimates of the expansion time in the BUU 
framework^ are surely larger than 100 fm/c. Hence the expansion process is expected 
to hardly affect the dynamics of fragment formation, while it could have some effects 
on the successive stages of the reaction. A realistic investigation of this following 
regime can be approached only by a fully dynamical simulation and it is beyond the 
scope of the present work. 

Now we discuss the reliability of the numerical estimates for A in our calculations. 
In fig. 9 we show the dependence of A on the initial distance do for average densities 
p/po=0.4, 0.5. There is a strong dependence of A on do for large values, which however 
saturates for values smaller than ~ 5. • 10~ 3 fm~ 2 In the calculations displayed in fig. 8 
we actually used a value do < 10~ 4 /m~ 2 ,as in ref.B, which lies in the saturation 
zone and therefore assures a stable numerical result. Last but not least, we also 
checked that A's do not change varying the size of the torus or changing the time step 
of the integration^. In fig. 10 the dependence of A on the initial nu is displayed for 



p/po=0.4,0.5. Even in this case the insensitivity to the variation of the wave number 
assures that the value of the Lyapunov exponent is numerically stable. It is important 
to notice that A is different from the growth rate calculated in the linear response 
theory, where one should expect a linear dependence with the largest mode having 
the fastest growth. Only if the system is linearly unstable the Lyapunov exponent 
and the growth rate do coincide. This is not true in our case where the latter is 
smaller by almost a factor of two. 

Up to now the momentum distribution adopted was always that one of a Fermi 
gas at temperature T=3 MeV. We varied also this temperature to see the effect on 
A. In fig. 11 the dependence of A on the temperature is shown. We found a small 
dependence on the temperature which indicates that thermal motion seems to reduce 
only slightly the degree of chaoticity. This is due to the shrinking of the spinodal 
zone at increasing temperatures, where vaporization processes start to compete with 
instabilities. 

The calculations discussed up to now refer only to the Vlasov equation. One 
could wonder how much the inclusion of two-body collisions can modify the discussed 
scenario. To this purpose, we solved the Vlasov-Nordheim equation numerically and 
calculated the Lyapunov exponents also in this case. In fig. 12 we show the comparison 
of X(t) calculated with and without two-body collisions. The figure shows that the 
inclusion of the collision integral (4) in the r.h.s of eq. (3) influences only slightly 
A reducing its value. This reduction is of the order of 30% for p/p =0.5, but only 
of 10% for p/p =0.4Ll. It is not difficult to understand this result by considering the 
diluteness of the fermion gas inside the spinodal region. 
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4.2. Time scales 

In this section we discuss the interplay of the different time scales occurring in 
fragment formation. 

In fig. 13 we show the comparison of the Lyapunov time T c h aos = 1/A and the time 
needed to form the primary fragments Tf rag , corresponding to the time end of the 
plateau. The former is always smaller than the latter. It is important to notice that 
Tf ra g is analogous to the trapping time T discussed in section 1. Hence according 
to the criterium there introduced for the simple example, the relation Tf rag > T c h aos 
implies that, during the interval of time in which fragments are formed, deterministic 
chaos can fully develop and characterizes the clusterization. 

An additional evidence of that can be also seen in fig. 7, where some modes, despite 
the overall growth of density fluctuation, can have amplitudes oscillating in time. This 
behaviour is strongly reminiscent of the erratic bouncing of the particles inside the 
potential of the example discussed in section 1. Of course the dynamics is followed in a 
finite time interval and therefore in our case chaoticity is only a transient phenomenon. 
Poincare maps like in ref.B3 cannot be drawn. We believe however that the evidence 
found are strong enough to justify at least a strong analogy of the above discussed 
dynamics with the more familiar low- dimensional chaoticity in closed and scattering 
systems H. 

It can be instructive to compare the characteristic times T c h aos = h/X, which 
defines the time scale of the divergence between mean field trajectories, with the 
single particle characteristic time r sp = h/Ep(p), being Ep the Fermi energy at the 
given density. This is done in table 1 for a set of densities and for a temperature T=3 
MeV. One can see that for densities p < 0.4p the divergence time is smaller than 
the single particle time. In other words the motion associated with the mean field is 
faster than that of the particles moving inside. Therefore in this region the notion 
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itself of mean field ceases to have validity, being the gas too dilute. 



5. Chaoticity implications for multifragmentation 

In this section we discuss what are the consequences of a chaotic dynamics in 
the spinodal zone for the final fragments formation. In our schematic description 
we are neglecting the initial compression induced by the heavy ion collision which 
likely leads the composite system well inside the spinodal zone. However one could 
simulate the uncertainty in the initial conditions by considering in our model a small 
and random initial perturbation in the density profile. This is also another way of 
studying the response of our system inside the spinodal region. In particular we start 
the time evolution giving to each cell along the x-direction a random shift around 
the average density whose strength is 1% of the average. Also with this initialization 
the time evolution is non-linear and show the same features already discussed in the 
previous section for a sinusoidal shape. This behaviour is shown in figs. 14 and 15 
where the time evolution of two random initialized events are compared by displaying 
the density profiles and the power spectra at different times. In fig. 14 the scale at 
t=0 is different to magnify the random initial shape. The power spectra of this shape 
is rather flat as fig. 15 illustrates. During the first evolution a wide range of modes 
are privileged due to the finite interaction range and then among these an erratic 
and hardly predictable mixing occurs. The dynamics of the single modes exhibits 
the same non-linear trend illustrated in fig. 7 a. We have checked that calculating the 
Lyapunov exponent with a random initialization one gets values identical to those 
shown previously. 

The initial random shape is on one hand more realistic, simulating the missing 
dynamics, and on the other hand it gives the opportunity to sample many initial 
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conditions. In this perspective we performed 108 runs for an average density p/po=0.5 
and extracted the fragment size distribution out of each one. In our case we do not 
have real fragments but a natural criterium is to consider those contiguous cells whose 
density exceeds a threshold value as forming a cluster. We adopted a freeze-out time 
which is slightly larger than that shown in fig. 13 and defined as Tf rag and is equal to 
120 fm/c. This was necessary because the dynamics is somehow slower when many 
modes are coupled since the beginning. Therefore the time for fragment formation 
shown in fig. 13 should be considered a minimum estimate. In fig. 16 we show the 
average fragment size distribution P(s) for three different threshold cuts equal to 
5,7 and 10 % of the normal density po- As one can see the three distributions are 
rather large and do not differ significantly between each other. In the same figure 
two fits for each distribution are shown. The dashed curve is an exponential fit, i.e. 
P(s) ~ e~ TS , while the full curve corresponds to a power law fit, i.e. P(s) ~ s~ T . The 
corresponding r and the relative errors are reported in the figure. The fits shown in 
the figure were performed up to a maximum value s max = 20. In table 2 we report the 
values obtained, together with the respective x 2 , by fitting the three distributions up 
to s max = 20, 30. In all cases the power law fits have a smaller \ 2 than the exponential 
ones. Though the difference is small, a power law fit always seems more appropriate. 
The value of the r power which one can extract from the power law fits oscillates around 
2. In a preliminary analysis of this kind we took into account 100 events considering 
a freeze-out time smaller and a larger cutH However a power law fit gave also in that 
case a similar value for r power . Thus the more refined present analysis does not change 
the previous results. 

One could be tempted to compare these calculations to those published recently 

no 
concerning experimental data and theoretical considerations Q'c3. In fact, the large 

distributions shown in fig. 16 remind those observed experimentally and the exponents 
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found are also very close to those of percolation at criticalitya Actually, these similar- 
ities could be accidental. Our model is at the moment too schematic and the criteria 
adopted are rather arbitrary to draw final conclusions. However, we think it is very 
important the fact that a deterministic chaotic behavior allows for final large distri- 
butions which resemble those found in experimental data. On the contrary we note 
that within a linear unstable evolutionE-M the most unstable mode has the largest 
growth rate and wins always over the others. Therefore, in the linear regime, the path 
followed by the system in its dynamical evolution from a random small initialization 
towards a macroscopic fragmentation is always the same, apart from a very small 
spreading width. Hence one obtains a strongly peaked and gaussian size distribution. 
Such a kind of multifragmentation has not yet been observed experimentally. 

In concluding this section we claim that a deterministic chaotic mechanism is 
the most natural explanation for the multifragmentation data. In a chaotic regime 
even selecting events with the same impact parameter, energy, temperature, etc. the 
small uncertainties which will always be present would lead inevitably to very large 
fluctuations with mass and charge distributions ranging from an exponential to a 
power law. Chaoticity does not exclude phase transitions or statistical hypotheses. 
On the other hand it seems likely to be a general feature underlying these scenarios. 



6. Conclusions 

The growth of density fluctuations, in the spinodal region of nuclear matter EOS, 
has been advocated as the mechanism of fragment formation in the multifragmenta- 
tion process observed in heavy ion collisions. According to this scenario, the nuclear 
fragment formation is similar to the process of droplet formation in an over-saturated 
classical vapour. We have studied the dynamics of such a process by solving the 
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Vlasov-Nordheim equation numerically on a lattice. The nuclear system in the spin- 
odal region was schematized by a two-dimensional fluid confined inside a fixed torus. 
The initially homogeneous system was perturbed by adding a small density fluctu- 
ation and the dynamics was followed up to the time of fragment formation, when 
several well separated density humps are apparent. Several characteristics, typical of 
chaotic dynamics, have been identified and studied. 

i) Strong sensitivity to the initial conditions is apparent in the evolution of the system. 
Small variations of the initial conditions lead to large deviations in the final pattern 
of the density profile and therefore of the fragment distribution. 
ii) The amplitudes of the modes increase, in general, non-exponent ially, with possible 
oscillating behaviour, at least for the wavelengths relevant for the fragment formation 
( A > 6fm ). This indicates that the modes are coupled and the dynamics is highly 
non-linear. 

Hi) The analysis in terms of Lyapunov exponents gives values different from zero 
and independent from the wavelength of the initial perturbation. This is in sharp 
contrast with the expectations for a linear dynamics, for which the inverse Lyapunov 
exponent coincides with the growth rate of the fluctuation and is proportional to the 
momentum of the mode. 

iv) The growth rate of the density fluctuations is longer than the inverse average 
Lyapunov exponent, and therefore chaotic dynamics has time to develop during the 
process of fragment formation. This criteria of identifying chaotic dynamics has been 
illustrated with a simple example in Sec. 1 

Points i-iv give strong evidences of a chaotic behaviour of the dynamics in the 
fragment formation process. Of course such a process takes place in a finite time 
interval, of the order of 100 fm/c, and therefore the traditional method of drawing 
Poincare maps to identify possible chaotic dynamics, cannot be used here. The sit- 
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uation, however, is similar to the case of chaotic scattering, and we believe that the 
evidences are strong enough to justify the conclusion of a chaotic dynamics of the 
process. 

This conclusion gives support to the statistical models used to analyze the mul- 
tifragmentation data, since chaotic dynamics entails a tendency of the system to fill 
uniformly the available phase space. In fact, if the dynamics is chaotic, strong fluctu- 
ations are expected from one collision event to another, each one ending in a different 
region of the available phase space. On the average, therefore, if chaoticity is strong 
enough, the population of the final channels will be dominated by phase space. 

Finally, a preliminary analysis of the fragment size distribution shows an approxi- 
mate power law. It has to be pointed out that the presence of chaotic dynamics does 
not exclude the possibility that the system passes through a phase transition. On the 
other hand a power law does not necessarily imply a phase transition. 

More detailed and realistic investigations have to be performed in order to con- 
firm this appealing scenario. It should be also clear that we have not presented a 
new model, but rather a new perspective in nuclear dynamics which seems more 
appropriate and promising. 
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TABLES 



P/Po 


A 


Tchaos 


T~sp 




c/fm 


fm/c 


fm/c 


0.7 


0. 


OC 


7.8 


0.6 


2.8 -10~ 2 


28.6 


9.1 


0.5 


6.8 • 10~ 2 


14.7 


10.9 


0.4 


9.6 • lO- 2 


10.4 


13.6 


0.3 


0.10 


10. 


18.1 


0.2 


0.10 


10. 


27.2 



Table 1 For a temperature T = 3 MeV, the Lyapunov exponent A, as calculated in the text, and 
the corresponding characteristic time T c haos = 1/A are reported. For comparison the single particle 
time r sp = h/ Ep(p) for various values of the average density outside and inside the spinodal region 
is shown. 
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Smax — ^U 










(P/P0)cut 


Tpower 


Xpower 


Texpon 


Xexpon 


0.05 


1.48±0.16 


0.087 


0.10±0.01 


0.12 


0.07 


2.17±0.49 


0.80 


0.15±0.04 


0.96 


0.10 


2.17±0.34 


0.61 


0.15±0.03 


0.86 



Smax — "3" 










(P/P0)cut 


Tpower 


Xpower 


Texpon 


Xexpon 


0.05 


0.96±0.13 


0.62 


0.05±0.01 


0.73 


0.07 


1.62±0.22 


1.770 


0.09±0.01 


2.13 


0.10 


2.10±0.15 


1.010 


0.12±0.01 


1.32 



Table 2 The values of the exponential and power law fits to the fragment distribution shown in 
fig. 16 are shown. The fits were performed up to a maximum s value s max = 20, 30. 
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FIGURES 

Fig. 1 Final value g/ of the coordinate q\ for a given set of initial conditions, defined by the total 
energy E = 3.9 , and the constraints ^p\ + qf = 3.7, q<i = 0. The initial value qi of the coordinate 
<7i is taken randomly in the corresponding allowed interval. 

Fig. 2 The free energy per particle F and the pressure P are respectively shown in part (a) and 
(b) . Both quantities are calculated as function of the density p/po and for different temperatures T. 

Fig. 3 Time evolution of the density profile (a) for an average density p//°o = 0.8. The ini- 
tial harmonic perturbation is characterized by a node number n^ — 4. The corresponding power 
spectrum is shown in (b). 

Fig. 4 The same as fig. 3 but for nk = 5. 

Fig. 5 Time evolution of the density profile for an average density p//°o = 0.4. The examples 
displayed above correspond to three different initial conditions, obtained by changing the node 
number rife = 3, 4, 5. 

Fig. 6 Power spectra corresponding to the density profiles shown in fig. 5. 

Fig. 7 Time evolution of the strength of the main modes excited in the profile of fig. 5 for n^ = 5 
is shown as open diamonds in steps of 5 fm/c. The linear fits over the first four points are drawn 
as full lines. The corresponding wavelengths are also reported. 

Fig. 8 The quantity X(t) is shown as a function of time for different values of the average density 
reported also in the figure. The dashed line correspond to the largest Lyapunov exponent A. The 
full curve is only to guide the eye. 
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Fig. 9 The value of A is shown as a function of the initial distance do for different values of the 
initial density p/po = 0.4, 0.5. See text for further details. 

Fig. 10 The value of A is shown as a function of the node number n^ characterizing the initial 
harmonic oscillation. The squares correspond to the initial average density p/po — 0.4, while the 
diamonds correspond to p/po = 0.5. The dashed lines represent the average value. 

Fig. 11 The value of A is shown as a function of the temperature T for an initial average density 
p/po = 0.5. 

Fig. 12 Time evolution of A(i) for an average initial density p/po = 0.5 respectively without 
(crosses) and with the collision integral I[f] (squares). 

Fig. 13 The Lyapunov time T c haos = 1/A (diamonds) is shown in comparison to the time re- 
quested by the system in order to form fragments Tf rag (squares) at the various average densities. 
The initial condition is an harmonic perturbation of the initial density profile. 

Fig. 14 Time evolution of two random initial density profiles. The average density is p/po = 0.5. 
Please note the magnification of the initial scale. 

Fig. 15 Power spectra corresponding to the density profiles displayed in fig. 14. Also in this case 
the initial scale is magnified. 
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Fig. 16 Fragment size distributions P(s) obtained by considering 108 events at an average density 
p/po = 0.5. P(s) is normalized to the total number of events considered. Each event is generated 
through a random initialization of the kind shown in fig. 14. The three panels correspond to the 
different density cuts chosen in order to define the fragments. In the figure are also drawn the 
power law fit (solid line) and the exponential one (dashed line). The maximum value for the fit was 
Smax == 20. The values of the fitted parameters with the relative errors are reported too. See text 
for further details. 
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